Predation interactions co-occurrence modelling

Packages

Code

# Data manipulation package
library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(psych) # Procedures for Psychological, Psychometric, and Personality Research

# Modelling packages
library(nlme) # Linear and Nonlinear Mixed Effects Models
library(betareg) # Beta Regression
library(glmmTMB) # Generalized Linear Mixed Models using Template Model Builder
library(AICcmodavg) # Model Selection and Multimodel Inference Based on (Q)AIC(c)
library(MuMIn) # Multi-Model Inference

#Modelling utility packages
library(broom) # Convert Statistical Objects into Tidy Tibbles
library(broom.mixed) # Tidying Methods for Mixed Models
library(insight) # Easy Access to Model Information for Various Model Objects
library(modelbased) # Estimation of Model-Based Predictions, Contrasts and Means
library(parameters) # Processing of Model Parameters
library(performance) # Assessment of Regression Models Performance
library(DHARMa) # Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models
library(ggeffects) # Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs

# Grapich packages
library(GGally) # Extension to 'ggplot2'
library(patchwork) # The Composer of Plots


# Table output packages
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax
library(DT)


library(showtext)
font_add_google("Roboto", "Roboto")
showtext_auto()
showtext_opts(dpi = 300)

Variables description

Table 1: List of variables used to modelling spatial co-occurrence patterns of predator prey interactions

Variable	Description	Units	Source
SIF	Spatial species interaction factor. SIF<1: segregation, SIF=1: independence, SIF>1= aggregation	index	Each study
Ov_coeff	Temporal overlap coefficient. Ov_coeff =1: complete activity overlap, Ov_Coeff =0: No activity overlap	index	Each study
CBMratio	Competitors body mass ratio: ln( dominant competitor body mass/ subordinate competitor body mass)	index	COMBINE data base (Soria et al. 2021)
Dmass_cat	Dominant competitor body mass category: <2 kg small, >2 kg and < 5kg Small-Medium, >5 kg and < 15kg Medium-Large, > 15 kg Large	category	COMBINE data base (Soria et al. 2021)
D_diet	Dominant competitor diet category: Carnivore if > 80% of vertebrate diet composition, Herbivore if > 80% of plant diet composition, Insectivore if >80% of invertebrate diet composition, Omnivore rest of the species	category	PHYLACINE trait data base (Faurby et al. 2020)
Diet_dist	Dominant and subordinate competitor category diets: Same if pair of competitors shares de diet category, different if not	category	PHYLACINE trait data base (Faurby et al. 2020)
Dist_tax	Dominant and subordinate competitor taxonomic family distance: Same if the pair of spaces share the taxonomic family, Different if not	category	COMBINE data base (Soria et al. 2021)
Phy_dist	Pairwise mitochondrial DNA phylogenetic distances	mitochondrial DNA phylogenetic distances	Hassanin et al. (2021)
Abs_lat	Absolute latitude of each study when reported	coordinates	Each study
Continent	Continent where the study was carried out	category	Each study
Label	Study labeling	category	Each study

Data

The data correspond to studies evaluating the spatial and temporal co-occurrence of mammals of the order Carnivora. In this database, mammals of the order Carnivora were considered as predators and mammals of other orders as potential prey.

Spatial data

Spatial co-occurrence data of predators and prey mammals were measured using the species interaction factor (SIF). This is a parameter derived from multi-species occupancy models (Waddle et al. 2010; Richmond, Hines, and Beissinger 2010).

Code

# Load the spatial co-ocurrence data base of predator-prey interactions
spatP_db <- read_delim("Data/Model data/spatP_db.csv", delim = ",")[,-1] %>%
  # Select the variables
  select(SIF, Pred_sp, Pred_family, Pred_mass, Prey_sp, Prey_family, Prey_mass, mass_ratio, P_hunt, Lon, Lat, Lat_abs, Avg_dist, Locality, Label, Pred_evidence) %>% 
  # filter species pair without predator-prey relation
  filter(Pred_evidence == "Yes")
DT::datatable(spatP_db)

The Table 2 contain the general description of the numeric variables of spatial data base. The table was constructed with psych package (Revelle 2022)

Code

Spatnum_summary <- spatP_db %>% 
  select_if(is.numeric) %>% 
  describe(. ,fast = T) 

kbl(Spatnum_summary, caption = "General description of the spatial co-occurrence data base for predator-prey interactions", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Table 2: General description of the spatial co-occurrence data base for predator-prey interactions.

	vars	n	mean	sd	min	max	range	se
SIF	1	25	1.081	0.823	0.473	4.656	4.183	0.165
Pred_mass	2	25	37320.851	63890.577	800.000	240500.000	239700.000	12778.115
Prey_mass	3	25	35982.358	46833.523	66.900	135000.000	134933.100	9366.705
mass_ratio	4	25	0.991	2.113	-1.777	5.774	7.551	0.423
Lon	5	23	19.578	91.034	-146.668	113.333	260.002	18.982
Lat	6	23	10.838	36.833	-41.234	63.824	105.057	7.680
Lat_abs	7	23	32.751	18.922	15.402	63.824	48.422	3.946
Avg_dist	8	21	1.214	1.393	0.500	4.000	3.500	0.304

We use the GGAlly package to explore the visual relationship between numerical variables (Schloerke et al. 2021). We used pearson’s correlation coefficient

Code

(SpatP_cor <- select_if(spatP_db, is.numeric) %>% # select numeric variables
  ggpairs(.,
          # Correlation coefficient upper part
        upper = list(continuous= wrap("cor", method= "pearson", 
                                      digits=2, corSize= 80)),
        lower = list( continuous= "smooth")) +
  theme_bw()+
  theme(text=element_text(family = "Roboto")))

Figure 1: Spatial predator-prey co-occurrence corplot

Temporal data

The temporal data corresponds to the co-occurrence of Carnivora mammals with their mammal prey measured by overlap coefficient of kernel activity curves (Ridout and Linkie 2009).

Code

tempP_db <- read_delim("Data/Model data/tempP_db.csv", delim=",")[,-1] %>%
  select(Ov_coeff, Pred_sp, Pred_family, Pred_mass, Prey_sp, Prey_family, Prey_mass, mass_ratio, P_hunt, Lon, Lat, Lat_abs, Samp_dur, Locality, Label,Pred_evidence) %>% 
  # filter species pair without predator-prey relation
  filter(Pred_evidence == "Yes")
DT::datatable(tempP_db )

The Table 3 contain the general description of the numeric variables of spatial data base.

Code

Tempnum_summary <- tempP_db %>% 
  select_if(is.numeric) %>% 
  describe(. ,fast = T) 
  
kbl(Tempnum_summary, caption = "General description of the temporal co-occurrence data base for predator-prey interactions", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Table 3: General description of the temporal co-occurrence data base for predator-prey interactions.

	vars	n	mean	sd	min	max	range	se
Ov_coeff	1	693	0.602	0.190	0.040	0.930	0.890	0.007
Pred_mass	2	685	44417.382	41252.700	110.330	240500.000	240389.670	1576.185
Prey_mass	3	693	71389.240	190227.527	36.000	3220000.000	3219964.000	7226.146
mass_ratio	4	685	0.734	1.669	-3.961	5.350	9.311	0.064
Lon	5	693	11.131	79.008	-123.084	149.183	272.267	3.001
Lat	6	693	16.199	21.837	-41.234	56.509	97.743	0.830
Lat_abs	7	693	22.753	14.873	1.686	56.509	54.823	0.565
Samp_dur	8	541	21.182	34.275	1.000	200.000	199.000	1.474

We use the GGAlly package to explore the visual relationship between numerical variables (Schloerke et al. 2021). We used pearson’s correlation coefficient

Code

(TempP_cor <- select_if(tempP_db, is.numeric)  %>% 
  ggpairs(.,
        upper = list(continuous= wrap("cor", method= "pearson", 
                                      digits=2, corSize= 80)),
        lower = list( continuous= "smooth")) +
  theme_bw()+
  theme(text=element_text(family = "Roboto")))

Figure 2: Temporal predator-prey co-occurrence correlogram

Modelling procedure

Standardize numeric covariates

The numerical covariates were standardized to mean 0 and standard deviation 1, which facilitates model fitting and direct comparison of regression coefficients. In addition, pairs of species that do not contain information on any of the variables (cells with NAs) were eliminated.

Because we aimed to evaluate whether the patterns found varied depending on the predator family, we made subsets of the databases. For spatial co-occurrence no subset were made. For temporal co-occurrence we subset the families Felidae and Canidae.

Code

scale_this <- function(data){scale(data) %>% as.vector() }

SP_vars <- c("mass_ratio", "Lat_abs", "Avg_dist")
TP_vars <- c("mass_ratio", "Lat_abs", "Samp_dur")

# All data Spatial
spatP_modsdf <- spatP_db %>%
  drop_na() %>% 
  mutate(across(all_of(SP_vars), scale_this)) 
dim(spatP_modsdf)

[1] 21 16

Code

# All data Temporal
tempP_modsdf <- tempP_db %>%
  drop_na() %>% 
  mutate(across(all_of(TP_vars), scale_this))
dim(tempP_modsdf)

[1] 538  16

Code

# Felidae Temporal
tempP_modsdf_F <- tempP_modsdf %>%
   filter(Pred_family == "Felidae")
dim(tempP_modsdf_F)

[1] 407  16

Code

# Canidae Temporal
tempP_modsdf_C <- tempP_modsdf %>%
   filter(Pred_family == "Canidae")
dim(tempP_modsdf_C)

[1] 85 16

Spatial modelling procedure

Influential observations

Identify the presence of outliers or observations that may influence the models. For the spatial co-occurrence data we fit a general linear model with all covariates and interactions. We then checked for the presence of extreme data using Laverage distance plots from performance package (Lüdecke et al. 2021).

Code

# Fit model with all data
SPmods_outl <- lm(SIF ~ mass_ratio+ I(mass_ratio^2)+ P_hunt+ Avg_dist, data= spatP_modsdf)

# Check outliers for each model
SP_out <- check_outliers(SPmods_outl)

# Leverage plot for all data model
plot(SP_out)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

Figure 3: Spatial co-occurrence leverage distance

Collinearity

We checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value > 10 indicates a high collinearity of the variables (Alain F. Zuur, Ieno, and Elphick 2010).

To calculate the VIF we use the Performance package (Lüdecke et al. 2021).

Code

spatP_coll <- lm(SIF ~ mass_ratio+ P_hunt+ Avg_dist, data= spatP_modsdf)

spatP_colltab <- check_collinearity(spatP_coll)


plot(spatP_colltab)+
  labs(subtitle = "All data Spatial model")+
  theme(legend.position = "none")

Figure 4: Spatial co-occurrence collinearity

Model assumptions

We first check the assumptions of the model. To do this we fit the most complex model containing all available variable interactions and random factors. We then visually evaluated the normality of the residuals, linearity and homogeneity of variance using the performance package (Lüdecke et al. 2021).

Code

Spat_diag_m <-lm(SIF ~ mass_ratio+ I(mass_ratio^2)+ P_hunt+  Avg_dist, 
                 random = ~1|Label/Locality,
                 method = "REML",
                 data= spatP_modsdf)

Spat_diagnostic <- check_model(Spat_diag_m, 
                               check = c("homogeneity", "linearity", "qq", "normality", "reqq"))
Spat_diagnostic 

No significant deviations from assumptions were detected.

Random structure

Following the Alain F. Zuur et al. (2009a) protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package (Mazerolle 2023).

Code

# Fit model without random effects
Spat_r0 <- glmmTMB(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+ Avg_dist,
                   data=spatP_modsdf, 
                   family=gaussian(link = "identity"), 
                   REML = T)

# Fit model with label as random effect
Spat_r1 <- glmmTMB(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+  Avg_dist +(1|Label),
                   data=spatP_modsdf, 
                   family=gaussian(link = "identity"), 
                   REML = T)

# Fit model with label and Locality as random effects
Spat_r2 <- glmmTMB(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+  Avg_dist +(1|Label/Locality),
                   data=spatP_modsdf, 
                   family=gaussian(link = "identity"), 
                   REML = T)


# Rank the models with the AIC
spat_AICtab <- aictab(cand.set = list(Spat_r0, Spat_r1, Spat_r2),
                       modnames = c("no random",
                                    "Label random", 
                                    "Label/Locality random"),
                                    second.ord = F)

kbl(spat_AICtab, caption = "Random-effects structure selection table for spatial co-occurrence data", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Table 4: Random-effects structure selection table for spatial co-occurrence data

Modnames	K	AIC	Delta_AIC	ModelLik	AICWt	LL	Cum.Wt
no random	7	19.478	0	1.000	0.665	-2.739	0.665
Label random	8	21.478	2	0.368	0.245	-2.739	0.910
Label/Locality random	9	23.478	4	0.135	0.090	-2.739	1.000

The Table 4 suggests that models without random effects is best supported. To visualize the variation of the random parameters of each group, we use the estimate_grouplevel function of the modelbased package (Makowski et al. 2020).

Code

estimate_grouplevel(Spat_r1) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

There is no significant variation from each group level in random effects.

Fixed predictor selections

To select the fixed effects structure we also used the Akaike information criterion (AIC). Models with a delta AIC <2 are considered equally plausible (Burnham and Anderson 2002). Because we are going to choose the fixed effects structure, we refit the previously identified model without using the restricted maximum likelihood (REML) (Alain F. Zuur et al. 2009b). because of the amount of data in the spatial co-occurrence model, we only consider in this case additive interactions of the variables. To do this we used the dredge function of the MuMIn package (Barton 2020).

All data

Code

SP_global <- lm(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+ Avg_dist,
                   data=spatP_modsdf, 
                   na.action = "na.fail")


SP_selec <- dredge(SP_global, rank = "AIC", m.lim= c(NA,3))

Code

kbl(SP_selec, 
    caption = "Fixed-effects structure selection table for spatial co-occurrence data", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Fixed-effects structure selection table for spatial co-occurrence data

	(Intercept)	Avg_dist	mass_ratio	I(mass_ratio^2)	P_hunt	df	logLik	AIC	delta	weight
15	0.733	NA	-0.174	0.082	+	6	8.791	-5.582	0.000	0.27396412
7	0.840	NA	-0.098	0.069	NA	4	6.013	-4.026	1.556	0.12582834
1	0.906	NA	NA	NA	NA	2	3.745	-3.490	2.092	0.09626930
11	0.809	NA	-0.095	NA	+	5	6.614	-3.227	2.355	0.08440608
3	0.906	NA	-0.053	NA	NA	3	4.466	-2.932	2.650	0.07282897
5	0.878	NA	NA	0.029	NA	3	4.073	-2.147	3.435	0.04918139
8	0.839	-0.015	-0.104	0.070	NA	5	6.069	-2.138	3.444	0.04896721
9	0.843	NA	NA	NA	+	4	5.025	-2.049	3.533	0.04683626
12	0.793	0.040	-0.082	NA	+	6	6.966	-1.931	3.651	0.04415665
10	0.809	0.063	NA	NA	+	5	5.853	-1.706	3.876	0.03945156
2	0.906	0.014	NA	NA	NA	3	3.792	-1.585	3.997	0.03713302
4	0.906	-0.007	-0.056	NA	NA	4	4.477	-0.954	4.628	0.02708790
6	0.876	0.018	NA	0.031	NA	4	4.157	-0.315	5.267	0.01967567
13	0.834	NA	NA	0.014	+	5	5.106	-0.213	5.369	0.01869998
14	0.803	0.062	NA	0.012	+	6	5.920	0.161	5.743	0.01551355

Confidence intervals explorations

Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models (Sutherland et al. 2023). Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference (Arnold 2010).

Code

# Function to get 85%IC
get_ci <- function(mod_name, mods){
  ci_df <- parameters::model_parameters(mods, ci = 0.85) %>% 
    select(Parameter, Coefficient, CI_low, CI_high, Component) %>% 
    filter(Component == "conditional") %>% 
    filter(Parameter != "(Intercept)") %>%
    mutate(Model = mod_name,
           Informative = case_when(
    CI_low < 0 & CI_high < 0 ~ "yes",
    CI_low > 0 & CI_high > 0 ~ "yes",
    TRUE ~ "no")) 
  return(ci_df)   
}


get_ci_lm <- function(mod_name, mods){
  ci_df <- parameters::model_parameters(mods, ci = 0.85) %>% 
    select(Parameter, Coefficient, CI_low, CI_high) %>% 
    filter(Parameter != "(Intercept)") %>%
    mutate(Model = mod_name,
           Informative = case_when(
    CI_low < 0 & CI_high < 0 ~ "yes",
    CI_low > 0 & CI_high > 0 ~ "yes",
    TRUE ~ "no")) 
  return(ci_df)   
}

# Function to plot

ci_plot <- function(ci_df){
ggplot(ci_df, aes(x=Coefficient, y= Parameter))+
  geom_pointrange(aes(xmin=CI_low, xmax= CI_high,
                      col= Model, linetype= Informative),
                  position = position_dodge2(0.5), linewidth= 1)+
  scale_linetype_manual(values = c("no"="dashed", "yes"= "solid"))+
  geom_vline(xintercept = 0, linetype= "dashed")+
  scale_color_viridis_d()+
  labs(caption = "estimates with 85% ci intervals",
       col= "Model ID")+
  theme_bw()+
    theme(text=element_text(family = "Roboto"))
  
}

Code

# get best models for all data 
SP_best_mods <- get.models(SP_selec, subset = delta <2)

# Apply the function to obtain the table of coefficients of selected the models
SP_best_ci <- map2_df(names(SP_best_mods), SP_best_mods, get_ci_lm) 

ci_plot(SP_best_ci)+ labs(title= "All data")

All data best model 85% CI

Temporal modelling procedure

Influential observations

Identify the presence of outliers or observations that may influence the models. For the temporal co-ocurrence data we fit a beta family generalized linear model with all covariates and interactions (Cribari-Neto and Zeileis 2010). We then checked for the presence of extreme observations using Cook’s distance plots from performance package (Lüdecke et al. 2021).

Code

# Models for each data
TPmods_outl <- betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, data= tempP_modsdf)

TPmods_outl_F <- betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, data= tempP_modsdf_F)

TPmods_outl_C <-betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, data= tempP_modsdf_C)


TP_out <- check_outliers(TPmods_outl)
TP_out_F <- check_outliers(TPmods_outl_F)
TP_out_C <- check_outliers(TPmods_outl_C)

# Cook's distance plots
plot(TP_out)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

plot(TP_out_F)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

plot(TP_out_C)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

Outliers of all data temporal model

Outliers of Felidae as predator Temporal model

Outliers of Canidae as predator Temporal model

Temporal co-occurrence cook distance distance

A possible outlier was detected for the Canidae data as predator. It corresponds to observation 27 which exceeded the 0.85 Cook’s distance threshold.

Collinearity

We checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value > 10 indicates a high collinearity of the variables (Alain F. Zuur, Ieno, and Elphick 2010).

To calculate the VIF we use the Performance package (Lüdecke et al. 2021).

All species

Code

# Fit models

tempP_coll <- betareg(Ov_coeff ~mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur, data= tempP_modsdf)
tempP_coll_F <- betareg(Ov_coeff ~mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur, data= tempP_modsdf_F)
tempP_coll_C <- betareg(Ov_coeff ~mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur, data= tempP_modsdf_C)

# Check collinearity
tempP_colltab <- check_collinearity(tempP_coll)
tempP_colltab_F <- check_collinearity(tempP_coll_F)
tempP_colltab_C <- check_collinearity(tempP_coll_C)

# Vif plots
plot(tempP_colltab)+
  labs(subtitle = "All data Temporal model")+
  theme_bw()+
  theme(legend.position = "none",
        text=element_text(family = "Roboto"))

plot(tempP_colltab_F)+
  labs(subtitle = "Felidae Temporal model")+
  theme_bw()+
  theme(legend.position = "none",
        text=element_text(family = "Roboto"))

plot(tempP_colltab_C)+
  labs(subtitle = "Canidae Temporal model")+
  theme_bw()+
  theme(legend.position = "none",
        text=element_text(family = "Roboto"))

(a) VIF of all data model

(b) VIF of Felidae model

(c) VIF of Canidae model

Figure 5: Outliers model detection temporal models

Model assumptions

We verify the model assumptions by visual inspection of the simulated residuals from the DHARMa package (Hartig 2021). To do so we fit the more complex model which includes fixed variables and their interactions, as well as random effects.

All data

Code

Temp_diag <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

Temp_diag_res <- simulateResiduals(Temp_diag)
plot(Temp_diag_res)

Since deviations from homogeneity of variance were detected, we evaluated which variables are involved in the model inadequacies.

Code

testCategorical(Temp_diag, tempP_modsdf$P_hunt)
testQuantiles(Temp_diag, tempP_modsdf$mass_ratio)
testQuantiles(Temp_diag, tempP_modsdf$Lat_abs)
testQuantiles(Temp_diag, tempP_modsdf$Samp_dur)

$uniformity
$uniformity$details
catPred: ambush

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.078041, p-value = 0.01688
alternative hypothesis: two-sided

------------------------------------------------------------ 
catPred: group

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.125, p-value = 0.4413
alternative hypothesis: two-sided

------------------------------------------------------------ 
catPred: pursuit

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.090531, p-value = 0.398
alternative hypothesis: two-sided


$uniformity$p.value
[1] 0.0168787 0.4413066 0.3979914

$uniformity$p.value.cor
[1] 0.0506361 0.7959827 0.7959827


$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)  
group   2  2.9148 0.05508 .
      535                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.01137
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.01306
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.0075
alternative hypothesis: both

Uniformity of variance test of Hunting strategy

Uniformity of variance test of mass ratio

Uniformity of variance test of absolute latitude

Uniformity of variance test of sampling duration

Code

Temp_diag2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf,
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), 
                   REML = T)

Temp_diag_res2 <- simulateResiduals(Temp_diag2)
plot(Temp_diag_res2)

Felidae

Code

Temp_diag_F <- glmmTMB(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), 
                   REML = T)


Temp_diag_res_F <- simulateResiduals(Temp_diag_F)
plot(Temp_diag_res_F)

Canidae

Code

Temp_diag_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality),
                   data=tempP_modsdf_C, 
                   family=beta_family(), 
                   REML = T)

Temp_diag_res_C <- simulateResiduals(Temp_diag_C)
plot(Temp_diag_res_C)

Code

Temp_diag_C2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality),
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), 
                   REML = T)

Temp_diag_res_C2 <- simulateResiduals(Temp_diag_C2)
plot(Temp_diag_res_C2)

Code

testCategorical(Temp_diag_C2, tempP_modsdf_C$P_hunt[-27])
testQuantiles(Temp_diag_C2, tempP_modsdf_C$mass_ratio[-27])
testQuantiles(Temp_diag_C2, tempP_modsdf_C$Lat_abs[-27])
testQuantiles(Temp_diag_C2, tempP_modsdf_C$Samp_dur[-27])

$uniformity
$uniformity$details
catPred: group

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.16372, p-value = 0.4185
alternative hypothesis: two-sided

------------------------------------------------------------ 
catPred: pursuit

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.098545, p-value = 0.6595
alternative hypothesis: two-sided


$uniformity$p.value
[1] 0.4185158 0.6594822

$uniformity$p.value.cor
[1] 0.8370315 0.8370315


$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.5443 0.4628
      82               

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.397
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.9715
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.5324
alternative hypothesis: both

Uniformity of variance test of Hunting strategy

Uniformity of variance test of mass ratio

Uniformity of variance test of absolute latitude

Uniformity of variance test of sampling duration

Random structure

Following the Alain F. Zuur et al. (2009a) protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the previously selected model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package (Mazerolle 2023).

Code

# All data
Temp_r0 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2,
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

Temp_r1 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ Lat_abs + (1|Label),
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

Temp_r2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality),
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

temp_AICtab <- aictab(cand.set = list(Temp_r0, Temp_r1, Temp_r2), 
                                  sort = T,
                                  modnames = c("no random",
                                               "Label random",
                                               "Label/Locality random")) 

Code

# Felidae
Temp_r0_F <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

Temp_r1_F <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ Lat_abs + (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

Temp_r2_F <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

temp_AICtab_F <- aictab(cand.set = list(Temp_r0_F, Temp_r1_F, Temp_r2_F), 
                                  sort = T,
                                  modnames = c("no random",
                                               "Label random",
                                               "Label/Locality random")) 

Code

# Canidae
Temp_r0_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, 
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), REML = T)

Temp_r1_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ Lat_abs + (1|Label), 
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), REML = T)

Temp_r2_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), REML = T)

temp_AICtab_C <- AICcmodavg::aictab(cand.set = list(Temp_r0_C, Temp_r1_C, Temp_r2_C), 
                                  sort = T,
                                  modnames = c("no random",
                                               "Label random",
                                               "Label/Locality random"))

Code

tempP_AICtab_all <- rbind(temp_AICtab, temp_AICtab_F, temp_AICtab_C)

kbl(tempP_AICtab_all, 
    caption = "Random-effects structure selection table for Temporal co-occurrence data", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F) %>%
  pack_rows("All data",1,3) %>%
  pack_rows("Felidae",4,6) %>% 
  pack_rows("Canidae",7, 9)

Random-effects structure selection table for Temporal co-occurrence data

	Modnames	K	AICc	Delta_AICc	ModelLik	AICcWt	LL	Cum.Wt
All data
2	Label random	26	-233.714	0.000	1.000	0.752	144.231	0.752
3	Label/Locality random	27	-231.497	2.217	0.330	0.248	144.231	1.000
1	no random	25	-216.630	17.085	0.000	0.000	134.584	1.000
Felidae
21	Label random	18	-252.375	0.000	1.000	0.750	145.069	0.750
31	Label/Locality random	19	-250.174	2.201	0.333	0.250	145.069	1.000
11	no random	17	-232.710	19.665	0.000	0.000	134.141	1.000
Canidae
12	no random	17	19.557	0.000	1.000	0.738	11.858	0.738
22	Label random	18	21.966	2.409	0.300	0.221	12.278	0.959
32	Label/Locality random	19	25.318	5.761	0.056	0.041	12.278	1.000

Code

estimate_grouplevel(Temp_r1) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

estimate_grouplevel(Temp_r1_F) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

estimate_grouplevel(Temp_r1_C) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

Group level estimates for all temporal data

Group level estimates for Felidae as predator temporal data

Group level estimates for Canidae as predator temporal data

Fixed predictor selections

To select the fixed effects structure we also used the Akaike information criterion (AIC). Models with a delta AIC <2 are considered equally plausible (Burnham and Anderson 2002). Because we are going to choose the fixed effects structure, we refit the previously identified model without using the restricted maximum likelihood (REML) (Alain F. Zuur et al. 2009b). Since co-occurrence patterns may be generated by the interaction of variables, we fit all possible combinations of variables, limiting to a maximum of three parameters per model. To do this we will used the dredge function of the MuMIn package (Barton 2020).

All data

Code

TP_global <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), 
                   REML = F,
                    na.action = "na.fail")

TP_selec <- dredge(TP_global, 
                   rank = "AIC", 
                   fixed = c("disp(mass_ratio)","disp(Lat_abs)","disp(Samp_dur)"), 
                   m.lim= c(NA,6))

Code

kbl(TP_selec, 
    caption = "Fixed-effects structure selection table for temporal co-occurrence data", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F)

Fixed-effects structure selection table for temporal co-occurrence data

	cond((Int))	disp((Int))	cond(Lat_abs)	cond(mass_ratio)	cond(I(mass_ratio^2))	cond(P_hunt)	cond(Samp_dur)	cond(Lat_abs:mass_ratio)	cond(Lat_abs:P_hunt)	cond(Lat_abs:Samp_dur)	cond(mass_ratio:P_hunt)	cond(mass_ratio:Samp_dur)	cond(I(mass_ratio^2):Lat_abs)	cond(I(mass_ratio^2):mass_ratio)	cond(I(mass_ratio^2):P_hunt)	cond(I(mass_ratio^2):Samp_dur)	cond(P_hunt:Samp_dur)	disp(Lat_abs)	disp(mass_ratio)	disp(Samp_dur)	df	logLik	AIC	delta	weight
74	0.451	1.911	0.085	NA	NA	+	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	0.026	-0.016	-0.071	11	169.361	-316.723	0.000	2.076382e-01
16409	0.461	1.921	NA	NA	NA	+	0.148	NA	NA	NA	NA	NA	NA	NA	NA	NA	+	0.027	-0.002	-0.036	11	169.228	-316.455	0.268	1.816312e-01
26	0.465	1.911	0.083	NA	NA	+	0.132	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.037	-0.006	-0.047	10	167.737	-315.474	1.248	1.112248e-01
25	0.470	1.910	NA	NA	NA	+	0.137	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.021	0.002	-0.060	9	166.710	-315.421	1.302	1.082803e-01
27	0.471	1.912	NA	-0.037	NA	+	0.139	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.026	-0.002	-0.059	10	167.205	-314.409	2.314	6.529908e-02
267	0.441	1.921	NA	-0.092	NA	+	NA	NA	NA	NA	+	NA	NA	NA	NA	NA	NA	-0.007	-0.015	-0.092	11	168.180	-314.359	2.364	6.368705e-02
10	0.442	1.908	0.088	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.026	-0.006	-0.063	9	165.971	-313.941	2.782	5.167866e-02
9	0.445	1.906	NA	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.006	0.002	-0.079	8	164.873	-313.746	2.977	4.686686e-02
29	0.477	1.911	NA	NA	-0.009	+	0.137	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.020	0.000	-0.059	10	166.769	-313.539	3.184	4.225497e-02
12	0.442	1.911	0.085	-0.032	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.029	-0.009	-0.063	10	166.347	-312.695	4.028	2.770819e-02
11	0.446	1.908	NA	-0.035	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.010	-0.001	-0.079	9	165.319	-312.637	4.085	2.692625e-02
14	0.450	1.909	0.088	NA	-0.009	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.024	-0.008	-0.061	10	166.037	-312.073	4.650	2.030532e-02
13	0.452	1.907	NA	NA	-0.009	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.004	0.000	-0.078	9	164.928	-311.857	4.866	1.822389e-02
15	0.454	1.909	NA	-0.036	-0.010	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.008	-0.005	-0.078	10	165.385	-310.771	5.952	1.058802e-02
4109	0.437	1.910	NA	NA	0.004	+	NA	NA	NA	NA	NA	NA	NA	NA	+	NA	NA	-0.002	-0.010	-0.084	11	166.135	-310.270	6.453	8.242532e-03
17	0.384	1.891	NA	NA	NA	NA	0.140	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.030	0.034	-0.065	7	160.425	-306.851	9.872	1.491318e-03
2055	0.370	1.910	NA	-0.135	-0.015	NA	NA	NA	NA	NA	NA	NA	NA	0.04	NA	NA	NA	0.036	0.025	-0.091	9	162.046	-306.092	10.631	1.020527e-03
21	0.403	1.894	NA	NA	-0.020	NA	0.141	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.026	0.027	-0.063	8	160.717	-305.435	11.288	7.347469e-04
18	0.375	1.892	0.043	NA	NA	NA	0.137	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.037	0.031	-0.061	8	160.710	-305.419	11.304	7.290206e-04
1	0.358	1.888	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.014	0.034	-0.086	6	158.613	-305.226	11.497	6.618472e-04
19	0.383	1.892	NA	-0.017	NA	NA	0.141	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.032	0.033	-0.065	8	160.537	-305.074	11.649	6.134844e-04
8213	0.399	1.896	NA	NA	-0.015	NA	0.167	NA	NA	NA	NA	NA	NA	NA	NA	-0.03	NA	0.026	0.031	-0.062	9	161.498	-304.997	11.726	5.902105e-04
146	0.367	1.891	0.043	NA	NA	NA	0.097	NA	NA	0.069	NA	NA	NA	NA	NA	NA	NA	0.038	0.032	-0.059	9	161.398	-304.797	11.926	5.340025e-04
22	0.394	1.894	0.045	NA	-0.021	NA	0.138	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.033	0.023	-0.059	9	161.023	-304.046	12.676	3.669839e-04
2	0.349	1.890	0.049	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.024	0.031	-0.080	7	158.954	-303.907	12.815	3.423555e-04
5	0.377	1.890	NA	NA	-0.020	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.011	0.027	-0.083	7	158.897	-303.795	12.928	3.236104e-04
1030	0.377	1.901	0.007	NA	-0.032	NA	NA	NA	NA	NA	NA	NA	0.044	NA	NA	NA	NA	0.028	0.022	-0.073	9	160.862	-303.724	12.999	3.123284e-04
23	0.402	1.895	NA	-0.018	-0.021	NA	0.142	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.028	0.025	-0.063	9	160.834	-303.667	13.056	3.035919e-04
20	0.374	1.893	0.042	-0.016	NA	NA	0.138	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.039	0.031	-0.061	9	160.800	-303.600	13.123	2.935664e-04
3	0.357	1.889	NA	-0.016	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.017	0.033	-0.086	7	158.707	-303.413	13.310	2.673875e-04
531	0.382	1.893	NA	-0.018	NA	NA	0.139	NA	NA	NA	NA	0.009	NA	NA	NA	NA	NA	0.032	0.032	-0.065	9	160.568	-303.137	13.586	2.328955e-04
6	0.368	1.893	0.051	NA	-0.021	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.020	0.023	-0.077	8	159.264	-302.527	14.196	1.717009e-04
4	0.348	1.891	0.047	-0.014	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.026	0.030	-0.080	8	159.027	-302.055	14.668	1.355573e-04
7	0.376	1.892	NA	-0.017	-0.020	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.013	0.025	-0.083	8	158.996	-301.992	14.730	1.314138e-04
36	0.342	1.891	0.054	-0.005	NA	NA	NA	-0.047	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.026	0.023	-0.078	9	159.903	-301.807	14.916	1.197700e-04
8	0.368	1.893	0.049	-0.015	-0.021	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.022	0.021	-0.077	9	159.342	-300.683	16.039	6.829379e-05

Felidae

Code

TP_Felidae <- glmmTMB(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label), 
                   data=tempP_modsdf_F,
                   family=beta_family(),
                   REML = F,
                   na.action = "na.fail")

TP_selec_F <- dredge(TP_Felidae, rank = "AIC", m.lim= c(NA,3))

Code

kbl(TP_selec_F, 
    caption = "Fixed-effects structure selection table for temporal co-occurrence data When Felidae is predator", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F)

Fixed-effects structure selection table for temporal co-occurrence data When Felidae is predator

	cond((Int))	disp((Int))	cond(Lat_abs)	cond(mass_ratio)	cond(I(mass_ratio^2))	cond(P_hunt)	cond(Samp_dur)	cond(Lat_abs:mass_ratio)	cond(Lat_abs:P_hunt)	cond(Lat_abs:Samp_dur)	cond(mass_ratio:P_hunt)	cond(mass_ratio:Samp_dur)	cond(I(mass_ratio^2):Lat_abs)	cond(I(mass_ratio^2):mass_ratio)	cond(I(mass_ratio^2):P_hunt)	cond(I(mass_ratio^2):Samp_dur)	cond(P_hunt:Samp_dur)	df	logLik	AIC	delta	weight
18	0.464	+	0.132	NA	NA	NA	0.105	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	155.852	-301.703	0.000	0.072760138
26	0.478	+	0.122	NA	NA	+	0.098	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	156.839	-301.677	0.026	0.071818780
10	0.466	+	0.127	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	155.749	-301.499	0.204	0.065698485
2	0.451	+	0.139	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	154.711	-301.422	0.281	0.063221822
20	0.464	+	0.123	-0.053	NA	NA	0.110	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	156.662	-301.323	0.380	0.060168047
12	0.466	+	0.119	-0.050	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	156.495	-300.991	0.713	0.050951490
4	0.450	+	0.131	-0.050	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	155.434	-300.869	0.835	0.047935630
267	0.461	+	NA	-0.076	NA	+	NA	NA	NA	NA	+	NA	NA	NA	NA	NA	NA	6	156.253	-300.506	1.197	0.039984181
8213	0.448	+	NA	NA	0.011	NA	0.150	NA	NA	NA	NA	NA	NA	NA	NA	-0.055	NA	6	156.203	-300.406	1.297	0.038040265
36	0.443	+	0.142	-0.060	NA	NA	NA	-0.058	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	156.172	-300.344	1.359	0.036875703
27	0.474	+	NA	-0.062	NA	+	0.108	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	155.997	-299.993	1.710	0.030943502
22	0.477	+	0.132	NA	-0.014	NA	0.104	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	155.941	-299.881	1.822	0.029257739
74	0.468	+	0.123	NA	NA	+	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	6	155.874	-299.747	1.956	0.027361915
146	0.466	+	0.132	NA	NA	NA	0.110	NA	NA	-0.012	NA	NA	NA	NA	NA	NA	NA	6	155.872	-299.744	1.959	0.027321172
25	0.473	+	NA	NA	NA	+	0.102	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	154.844	-299.688	2.016	0.026559184
19	0.460	+	NA	-0.062	NA	NA	0.118	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	154.843	-299.685	2.018	0.026527940
14	0.478	+	0.127	NA	-0.013	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	155.830	-299.660	2.043	0.026196281
6	0.465	+	0.139	NA	-0.016	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	154.829	-299.657	2.046	0.026161193
17	0.460	+	NA	NA	NA	NA	0.112	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	153.731	-299.461	2.242	0.023719992
11	0.461	+	NA	-0.059	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	154.701	-299.403	2.300	0.023034874
9	0.461	+	NA	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	153.671	-299.341	2.362	0.022338752
8	0.463	+	0.132	-0.049	-0.014	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	155.527	-299.055	2.648	0.019354848
3	0.445	+	NA	-0.059	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	153.449	-298.897	2.806	0.017890622
1	0.446	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	3	152.445	-298.891	2.812	0.017831724
1030	0.457	+	0.108	NA	-0.007	NA	NA	NA	NA	NA	NA	NA	0.034	NA	NA	NA	NA	6	155.232	-298.464	3.239	0.014402868
16409	0.476	+	NA	NA	NA	+	0.101	NA	NA	NA	NA	NA	NA	NA	NA	NA	+	6	155.102	-298.204	3.499	0.012649437
29	0.485	+	NA	NA	-0.013	+	0.102	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	154.916	-297.833	3.871	0.010505669
23	0.471	+	NA	-0.060	-0.013	NA	0.117	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	154.914	-297.829	3.874	0.010485735
531	0.459	+	NA	-0.063	NA	NA	0.114	NA	NA	NA	NA	0.014	NA	NA	NA	NA	NA	6	154.910	-297.821	3.883	0.010442640
21	0.473	+	NA	NA	-0.015	NA	0.112	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	153.837	-297.674	4.030	0.009702304
15	0.472	+	NA	-0.058	-0.012	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	6	154.762	-297.523	4.180	0.008999102
13	0.474	+	NA	NA	-0.014	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	153.761	-297.522	4.181	0.008992944
5	0.461	+	NA	NA	-0.017	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	152.578	-297.157	4.546	0.007492844
7	0.459	+	NA	-0.058	-0.015	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	153.547	-297.094	4.609	0.007263183
4109	0.477	+	NA	NA	-0.019	+	NA	NA	NA	NA	NA	NA	NA	NA	+	NA	NA	6	154.047	-296.094	5.609	0.004404904
2055	0.458	+	NA	-0.048	-0.014	NA	NA	NA	NA	NA	NA	NA	NA	-0.004	NA	NA	NA	6	153.559	-295.118	6.585	0.002704093

Canidae

Code

tempP_modsdf_C2 <- slice(tempP_modsdf_C, -27)

TP_Canidae <- betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+P_hunt + Lat_abs+Samp_dur)^2,
                         data=tempP_modsdf_C2, 
                         na.action = "na.fail")

TP_selec_C <- dredge(TP_Canidae, rank = "AIC", m.lim= c(NA,3))

Code

kbl(TP_selec_C, 
    caption = "Fixed-effects structure selection table for temporal co-occurrence data When Canidae is predator", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F)

Fixed-effects structure selection table for temporal co-occurrence data When Canidae is predator

	(Intercept)	Lat_abs	mass_ratio	I(mass_ratio^2)	P_hunt	Samp_dur	Lat_abs:mass_ratio	Lat_abs:P_hunt	Lat_abs:Samp_dur	mass_ratio:P_hunt	mass_ratio:Samp_dur	I(mass_ratio^2):Lat_abs	I(mass_ratio^2):mass_ratio	I(mass_ratio^2):P_hunt	I(mass_ratio^2):Samp_dur	P_hunt:Samp_dur	df	logLik	AIC	delta	weight
74	-0.336	0.539	NA	NA	+	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	5	20.551	-31.102	0.000	0.391806844
36	0.191	0.100	0.221	NA	NA	NA	-0.231	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	20.088	-30.176	0.926	0.246623687
1030	0.323	-0.071	NA	-0.142	NA	NA	NA	NA	NA	NA	NA	0.12	NA	NA	NA	NA	5	19.023	-28.047	3.055	0.085067398
26	-0.101	0.180	NA	NA	+	0.226	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	17.832	-25.663	5.438	0.025833351
18	0.098	0.176	NA	NA	NA	0.203	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	16.510	-25.019	6.082	0.018721043
22	0.178	0.187	NA	-0.053	NA	0.209	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	17.334	-24.667	6.435	0.015697255
16409	0.016	NA	NA	NA	+	0.483	NA	NA	NA	NA	NA	NA	NA	NA	NA	+	5	17.224	-24.447	6.654	0.014064257
25	0.016	NA	NA	NA	+	0.184	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	16.206	-24.413	6.689	0.013822389
2055	0.263	NA	-0.232	-0.045	NA	NA	NA	NA	NA	NA	NA	NA	0.061	NA	NA	NA	5	17.138	-24.276	6.826	0.012909008
1	0.186	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	2	14.044	-24.088	7.013	0.011753387
17	0.207	NA	NA	NA	NA	0.157	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	3	14.997	-23.994	7.107	0.011213252
20	0.123	0.163	0.072	NA	NA	0.227	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	16.970	-23.940	7.162	0.010911143
9	0.015	NA	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	3	14.960	-23.921	7.181	0.010808779
2	0.100	0.133	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	3	14.960	-23.920	7.182	0.010803512
10	-0.071	0.135	NA	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	15.905	-23.810	7.291	0.010228099
146	0.128	0.162	NA	NA	NA	0.200	NA	NA	0.075	NA	NA	NA	NA	NA	NA	NA	5	16.841	-23.682	7.420	0.009592431
19	0.227	NA	0.087	NA	NA	0.193	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	15.669	-23.339	7.763	0.008079372
6	0.175	0.143	NA	-0.050	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	15.668	-23.336	7.766	0.008068637
5	0.258	NA	NA	-0.045	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	3	14.600	-23.199	7.902	0.007534955
21	0.282	NA	NA	-0.046	NA	0.160	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	15.594	-23.188	7.914	0.007492804
29	0.092	NA	NA	-0.036	+	0.184	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	16.569	-23.138	7.963	0.007309855
14	0.013	0.142	NA	-0.040	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	16.374	-22.748	8.353	0.006013944
3	0.196	NA	0.056	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	3	14.331	-22.662	8.439	0.005761653
13	0.091	NA	NA	-0.036	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	15.318	-22.637	8.465	0.005688479
27	0.053	NA	0.033	NA	+	0.193	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	16.276	-22.552	8.550	0.005451713
4109	0.227	NA	NA	-0.099	+	NA	NA	NA	NA	NA	NA	NA	NA	+	NA	NA	5	16.264	-22.529	8.573	0.005388523
23	0.289	NA	0.077	-0.039	NA	0.192	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	16.110	-22.220	8.881	0.004619106
4	0.113	0.124	0.040	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	15.103	-22.206	8.895	0.004586076
11	0.016	NA	0.001	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	14.961	-21.921	9.181	0.003976818
12	-0.105	0.141	-0.026	NA	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	15.948	-21.896	9.205	0.003927985
7	0.260	NA	0.045	-0.040	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	4	14.789	-21.578	9.523	0.003350397
8	0.180	0.137	0.025	-0.047	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	15.724	-21.448	9.654	0.003138977
531	0.228	NA	0.088	NA	NA	0.194	NA	NA	NA	NA	0.002	NA	NA	NA	NA	NA	5	15.669	-21.339	9.763	0.002972685
8213	0.292	NA	NA	-0.052	NA	0.196	NA	NA	NA	NA	NA	NA	NA	NA	-0.023	NA	5	15.631	-21.263	9.839	0.002861382
15	0.087	NA	-0.004	-0.036	+	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	5	15.319	-20.639	10.463	0.002094630
267	0.125	NA	0.097	NA	+	NA	NA	NA	NA	+	NA	NA	NA	NA	NA	NA	5	15.182	-20.364	10.737	0.001826176

Confidence intervals explorations

Code

# Function to get 85%IC
get_cibeta <- function(mod_name, mods){
  ci_df <- parameters::model_parameters(mods, ci = 0.85) %>% 
    select(Parameter, Coefficient, CI_low, CI_high) %>% 
    filter(Parameter != "(Intercept)") %>%
    mutate(Model = mod_name,
           Informative = case_when(
    CI_low < 0 & CI_high < 0 ~ "yes",
    CI_low > 0 & CI_high > 0 ~ "yes",
    TRUE ~ "no")) 
  return(ci_df)   
}

Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models (Sutherland et al. 2023). Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference (Arnold 2010).

Code

# get best models for all data 
TP_best_mods <- get.models(TP_selec, subset = delta <2, REML = T)
# get best models for felidae subset data
TP_best_mods_F <- get.models(TP_selec_F, subset = delta <2, REML = T)
# get best models for Canidae subset data
TP_best_mods_C <- get.models(TP_selec_C, subset = delta <2)


# Apply the function to obtain the table of coefficients of selected the models
TP_best_ci <- map2_df(names(TP_best_mods), TP_best_mods, get_ci) 
TP_best_ci_F <- map2_df(names(TP_best_mods_F), TP_best_mods_F, get_ci)
TP_best_ci_C <- map2_df(names(TP_best_mods_C), TP_best_mods_C, get_cibeta)


ci_plot(TP_best_ci)+ labs(title= "All data")
ci_plot(TP_best_ci_F)+ labs(title= "Felidae as predator")
ci_plot(TP_best_ci_C)+ labs(title= "Canidae as predator")

All data best model 85% CI

Felidae data best model 85% CI

Canidae data best model 85% CI

Summary selected models

Note that some variables have coefficients whose 85% confidence intervals overlap 0. These variables are retained when they have interactions with other variables and these interactions do not overlap with 0.

Spatial co-occurrence models

Code

SP_mod1 <- lm(formula = SIF~ mass_ratio +I(mass_ratio^2)+  P_hunt,
                   data=spatP_modsdf)

SP_mod2 <- lm(formula = SIF~ mass_ratio +I(mass_ratio^2),
                   data=spatP_modsdf)


SP_mods <- list(SP_mod1, SP_mod2) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("SIF~ mass_ratio +I(mass_ratio)^2+P_hunt", 5),
                  rep("SIF~ mass_ratio +I(mass_ratio)^2", 3)),
         Component= "conditional",
         Effects= "fixed") %>% 
  select(-SE,-t,-df_error)

kbl(SP_mods, 
    caption = "Spatial co-occurrence selected models", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Spatial co-occurrence selected models

Parameter	Coefficient	CI	CI_low	CI_high	p	Model	Component	Effects
(Intercept)	0.733	0.85	0.624	0.843	0.000	SIF~ mass_ratio +I(mass_ratio)^2+P_hunt	conditional	fixed
mass_ratio	-0.174	0.85	-0.276	-0.073	0.020	SIF~ mass_ratio +I(mass_ratio)^2+P_hunt	conditional	fixed
I(mass_ratio^2)	0.082	0.85	0.018	0.147	0.073	SIF~ mass_ratio +I(mass_ratio)^2+P_hunt	conditional	fixed
P_huntgroup	-0.014	0.85	-0.267	0.239	0.935	SIF~ mass_ratio +I(mass_ratio)^2+P_hunt	conditional	fixed
P_huntpursuit	0.221	0.85	0.069	0.374	0.043	SIF~ mass_ratio +I(mass_ratio)^2+P_hunt	conditional	fixed
(Intercept)	0.840	0.85	0.753	0.927	0.000	SIF~ mass_ratio +I(mass_ratio)^2	conditional	fixed
mass_ratio	-0.098	0.85	-0.174	-0.021	0.072	SIF~ mass_ratio +I(mass_ratio)^2	conditional	fixed
I(mass_ratio^2)	0.069	0.85	0.008	0.130	0.108	SIF~ mass_ratio +I(mass_ratio)^2	conditional	fixed

Goftest of spatial selected models

Code

SP_best_goft_list <- list(SP_mod1, SP_mod2) 
SP_best_goft <- lapply(SP_best_goft_list, check_model, check = c("homogeneity", "linearity", "qq", "normality", "reqq"))

Temporal co-occurrence models

Code

TP_mod1 <-  glmmTMB(formula = Ov_coeff ~  Lat_abs *P_hunt + (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~ mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), REML = T)

TP_mod2 <-  glmmTMB(formula = Ov_coeff ~ P_hunt* Samp_dur + (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~ mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), REML = T)

TP_mod3 <-  glmmTMB(formula = Ov_coeff ~ P_hunt+ Samp_dur + (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~ mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), REML = T)

TP_mods <- list(TP_mod1, TP_mod2, TP_mod3) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("Ov_coeff ~  Lat_abs *P_hunt + (1|Label)", 11),
                  rep("Ov_coeff ~ P_hunt* Samp_dur + (1|Label)", 11),
                  rep("Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)", 9))) %>% 
  select(-SE, -z, -df_error, -Group)

Code

TP_mod1_F <-  glmmTMB(formula = Ov_coeff ~  Lat_abs + (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

TP_mod2_F <-  glmmTMB(formula = Ov_coeff ~  Lat_abs + Samp_dur+ (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)  

TP_mod3_F <-  glmmTMB(formula = Ov_coeff ~ mass_ratio *P_hunt + (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T) 

TP_mod4_F <-  glmmTMB(formula = Ov_coeff ~  I(mass_ratio^2)* Samp_dur+(1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

TP_mod5_F <-  glmmTMB(formula = Ov_coeff ~  mass_ratio+ Samp_dur +(1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

TP_mods_F <- list(TP_mod1_F, TP_mod2_F, TP_mod3_F, TP_mod4_F, TP_mod5_F) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("Ov_coeff ~  Lat_abs + (1|Label)", 4),
                  rep("Ov_coeff ~  Lat_abs + Samp_dur+ (1|Label)", 5),
                  rep("Ov_coeff ~ mass_ratio * P_hunt + (1|Label)", 6),
                  rep("Ov_coeff ~  I(mass_ratio^2) * Samp_dur+(1|Label)", 6),
                  rep("Ov_coeff ~  mass_ratio+ Samp_dur +(1|Label)", 5))) %>% 
  select(-SE, -z, -df_error,-Group)  

Code

TP_mod1_C  <-betareg(Ov_coeff ~ P_hunt * Lat_abs,
                         data=tempP_modsdf_C2)

TP_mod2_C  <-betareg(Ov_coeff ~ mass_ratio * Lat_abs,
                         data=tempP_modsdf_C2)

TP_mods_C <- list(TP_mod1_C, TP_mod2_C) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("Ov_coeff ~ P_hunt * Lat_abs", 4),
                  rep("Ov_coeff ~ mass_ratio * Lat_abs", 4)),
         Component= "conditional",
         Effects= "fixed") %>% 
  select(-SE, -z, -df_error)

Code

TP_mods_tab <- rbind(TP_mods, TP_mods_F, TP_mods_C)

kbl(TP_mods_tab, 
    caption = "Temporal co-occurrence selected models", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F) %>% 
    pack_rows("Temporal co-occurrence models -all data",1, nrow(TP_mods)) %>%
    pack_rows("Temporal co-occurrence models -Felidae as predator",1+ nrow(TP_mods),nrow(TP_mods)+ nrow(TP_mods_F)) %>%
  pack_rows("Temporal co-occurrence models -Canidae as predator",1+nrow(TP_mods)+ nrow(TP_mods_F), nrow(TP_mods)+ nrow(TP_mods_F)+ nrow(TP_mods_C))

Temporal co-occurrence selected models

Parameter	Coefficient	CI	CI_low	CI_high	p	Effects	Component	Model
Temporal co-occurrence models -all data
(Intercept)	0.444	0.85	0.354	0.535	0.000	fixed	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs	0.083	0.85	-0.012	0.178	0.209	fixed	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
P_huntgroup	-0.564	0.85	-0.769	-0.360	0.000	fixed	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
P_huntpursuit	-0.148	0.85	-0.318	0.022	0.209	fixed	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs:P_huntgroup	0.249	0.85	0.034	0.463	0.095	fixed	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs:P_huntpursuit	-0.146	0.85	-0.293	0.001	0.153	fixed	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
(Intercept)	1.911	0.85	1.822	2.000	0.000	fixed	dispersion	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
mass_ratio	-0.016	0.85	-0.108	0.075	0.798	fixed	dispersion	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Samp_dur	-0.072	0.85	-0.166	0.021	0.264	fixed	dispersion	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs	0.028	0.85	-0.075	0.130	0.695	fixed	dispersion	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
SD (Intercept)	0.313	0.85	0.238	0.411	NA	random	conditional	Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
(Intercept)	0.453	0.85	0.358	0.548	0.000	fixed	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntgroup	-0.438	0.85	-0.642	-0.235	0.002	fixed	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntpursuit	-0.154	0.85	-0.309	0.001	0.153	fixed	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
Samp_dur	0.147	0.85	0.037	0.258	0.054	fixed	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntgroup:Samp_dur	0.202	0.85	-0.101	0.505	0.338	fixed	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntpursuit:Samp_dur	-0.162	0.85	-0.294	-0.030	0.077	fixed	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
(Intercept)	1.919	0.85	1.830	2.009	0.000	fixed	dispersion	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
mass_ratio	-0.002	0.85	-0.092	0.089	0.975	fixed	dispersion	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
Samp_dur	-0.045	0.85	-0.134	0.045	0.473	fixed	dispersion	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
Lat_abs	0.024	0.85	-0.077	0.125	0.731	fixed	dispersion	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
SD (Intercept)	0.333	0.85	0.257	0.432	NA	random	conditional	Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
(Intercept)	0.463	0.85	0.368	0.557	0.000	fixed	conditional	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
P_huntgroup	-0.468	0.85	-0.667	-0.269	0.001	fixed	conditional	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
P_huntpursuit	-0.169	0.85	-0.324	-0.014	0.116	fixed	conditional	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Samp_dur	0.135	0.85	0.030	0.240	0.064	fixed	conditional	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
(Intercept)	1.910	0.85	1.822	1.999	0.000	fixed	dispersion	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
mass_ratio	0.001	0.85	-0.089	0.091	0.988	fixed	dispersion	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Samp_dur	-0.062	0.85	-0.152	0.027	0.316	fixed	dispersion	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Lat_abs	0.022	0.85	-0.077	0.122	0.745	fixed	dispersion	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
SD (Intercept)	0.329	0.85	0.256	0.423	NA	random	conditional	Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Temporal co-occurrence models -Felidae as predator
(Intercept)	0.444	0.85	0.354	0.535	0.000	fixed	conditional	Ov_coeff ~ Lat_abs + (1|Label)
Lat_abs	0.137	0.85	0.042	0.231	0.038	fixed	conditional	Ov_coeff ~ Lat_abs + (1|Label)
(Intercept)	7.977	0.85	7.205	8.831	NA	fixed	dispersion	Ov_coeff ~ Lat_abs + (1|Label)
SD (Intercept)	0.320	0.85	0.245	0.417	NA	random	conditional	Ov_coeff ~ Lat_abs + (1|Label)
(Intercept)	0.458	0.85	0.368	0.548	0.000	fixed	conditional	Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
Lat_abs	0.129	0.85	0.036	0.223	0.046	fixed	conditional	Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
Samp_dur	0.101	0.85	-0.001	0.204	0.154	fixed	conditional	Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
(Intercept)	7.973	0.85	7.202	8.826	NA	fixed	dispersion	Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
SD (Intercept)	0.309	0.85	0.235	0.408	NA	random	conditional	Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
(Intercept)	0.454	0.85	0.365	0.544	0.000	fixed	conditional	Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
mass_ratio	-0.075	0.85	-0.136	-0.015	0.073	fixed	conditional	Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
P_huntgroup	-0.244	0.85	-0.591	0.102	0.310	fixed	conditional	Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
mass_ratio:P_huntgroup	0.316	0.85	0.056	0.576	0.080	fixed	conditional	Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
(Intercept)	7.966	0.85	7.188	8.827	NA	fixed	dispersion	Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
SD (Intercept)	0.307	0.85	0.230	0.410	NA	random	conditional	Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
(Intercept)	0.443	0.85	0.342	0.544	0.000	fixed	conditional	Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
I(mass_ratio^2)	0.010	0.85	-0.040	0.060	0.777	fixed	conditional	Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
Samp_dur	0.147	0.85	0.042	0.253	0.045	fixed	conditional	Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
I(mass_ratio^2):Samp_dur	-0.053	0.85	-0.089	-0.018	0.032	fixed	conditional	Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
(Intercept)	7.997	0.85	7.220	8.858	NA	fixed	dispersion	Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
SD (Intercept)	0.316	0.85	0.240	0.416	NA	random	conditional	Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
(Intercept)	0.454	0.85	0.361	0.546	0.000	fixed	conditional	Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
mass_ratio	-0.061	0.85	-0.120	-0.002	0.139	fixed	conditional	Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
Samp_dur	0.115	0.85	0.012	0.219	0.108	fixed	conditional	Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
(Intercept)	7.990	0.85	7.215	8.847	NA	fixed	dispersion	Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
SD (Intercept)	0.326	0.85	0.250	0.425	NA	random	conditional	Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
Temporal co-occurrence models -Canidae as predator
(Intercept)	-0.336	0.85	-0.596	-0.076	0.063	fixed	conditional	Ov_coeff ~ P_hunt * Lat_abs
P_huntpursuit	0.675	0.85	0.351	0.999	0.003	fixed	conditional	Ov_coeff ~ P_hunt * Lat_abs
Lat_abs	0.539	0.85	0.305	0.774	0.001	fixed	conditional	Ov_coeff ~ P_hunt * Lat_abs
P_huntpursuit:Lat_abs	-0.620	0.85	-0.909	-0.331	0.002	fixed	conditional	Ov_coeff ~ P_hunt * Lat_abs
(Intercept)	0.191	0.85	0.029	0.353	0.089	fixed	conditional	Ov_coeff ~ mass_ratio * Lat_abs
mass_ratio	0.221	0.85	0.088	0.354	0.017	fixed	conditional	Ov_coeff ~ mass_ratio * Lat_abs
Lat_abs	0.100	0.85	-0.040	0.240	0.305	fixed	conditional	Ov_coeff ~ mass_ratio * Lat_abs
mass_ratio:Lat_abs	-0.231	0.85	-0.337	-0.124	0.002	fixed	conditional	Ov_coeff ~ mass_ratio * Lat_abs

Goftest of temporal selected models

Code

TP_best_goft_list <- list(TP_mod1, TP_mod2, TP_mod3) 
TP_best_goft <- lapply(TP_best_goft_list, simulateResiduals, plot=T)

Ov_coeff ~ Lat_abs * P_hunt + (1 | Label) goftest -all data

Ov_coeff ~ P_hunt * Samp_dur + (1 | Label) goftest - -all data

Ov_coeff ~ P_hunt + Samp_dur + (1 | Label) goftest- -all data

Code

TP_best_goft_list_F <- list(TP_mod1_F, TP_mod2_F, TP_mod3_F, TP_mod4_F, TP_mod5_F) 
TP_best_goft_F <- lapply(TP_best_goft_list_F, simulateResiduals, plot=T)

Ov_coeff ~ Lat_abs + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ Lat_abs + Samp_dur + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ mass_ratio * P_hunt + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ I(mass_ratio^2) * Samp_dur + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ mass_ratio + Samp_dur + (1 | Label) goftest -Felidae as predator

Code

TP_best_goft_list_C <- list(TP_mod1_C, TP_mod2_C) 
TP_best_goft_C <- lapply(TP_best_goft_list_C, plot)

Predictions

Code

# Spatial data

SP_pred_mass <- ggemmeans(SP_mod1, terms = c("mass_ratio[all]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(spatP_db$mass_ratio), "scaled:scale"))+ attr(scale(spatP_db$mass_ratio), "scaled:center"))

SP_pred_hunt <- ggemmeans(SP_mod1, terms = c("P_hunt"),
                          ci_level=.85,
                          ci.lvl = .85)

# Temporal data

TP_pred_lat1 <- ggemmeans(TP_mod1, 
                          terms = c("Lat_abs[-1.10:2.35 by=.2]",
                                    "P_hunt[pursuit]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lat2 <- ggemmeans(TP_mod1, 
                          terms = c("Lat_abs[-1.3:2.36 by=.2]",
                                    "P_hunt[ambush]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lat3 <- ggemmeans(TP_mod1, 
                          terms = c("Lat_abs[-1.15:2.37 by=.2]",
                                    "P_hunt[group]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lat <- rbind(TP_pred_lat1, TP_pred_lat2, TP_pred_lat3)

# PhuntxSampling
TP_pred_sxph1 <- ggemmeans(TP_mod2, 
                          terms = c("Samp_dur[-0.58:2.57 by=.2]",
                                    "P_hunt[pursuit]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

TP_pred_sxph2 <- ggemmeans(TP_mod2, 
                          terms = c("Samp_dur[-0.58:5.20 by=.2]",
                                    "P_hunt[ambush]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

TP_pred_sxph3 <- ggemmeans(TP_mod2, 
                          terms = c("Samp_dur[-0.53:2.14 by=.2]",
                                    "P_hunt[group]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

TP_pred_sxph <- rbind(TP_pred_sxph1, TP_pred_sxph2, TP_pred_sxph3)

# Felidae temporal data
#laitud
TP_pred_lat_F <- ggemmeans(TP_mod1_F, terms = c("Lat_abs"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

# masa * P_hunt

TP_pred_mxhunt_F1 <- ggemmeans(TP_mod3_F, 
                               terms = c("mass_ratio[-2.41:2.56 by=.2]", 
                                         "P_hunt[ambush]"),
                               ci_level=.85,
                               ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(tempP_db$mass_ratio), "scaled:scale"))+ attr(scale(tempP_db$mass_ratio), "scaled:center"))

TP_pred_mxhunt_F2 <- ggemmeans(TP_mod3_F, 
                               terms = c("mass_ratio[-1.48:1.57 by=.2]", 
                                         "P_hunt[group]"),
                               ci_level=.85,
                               ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(tempP_db$mass_ratio), "scaled:scale"))+ attr(scale(tempP_db$mass_ratio), "scaled:center"))


TP_pred_mxhunt_F <- rbind(TP_pred_mxhunt_F1, TP_pred_mxhunt_F2)

# Sampdur

TP_pred_samp_F <- ggemmeans(TP_mod5_F, 
                               terms = c("Samp_dur[-0.58:5.20 by=0.2]"),
                               ci_level=.85,
                               ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

# Canidae temporal data

TP_pred_lxh_C1 <-  ggemmeans(TP_mod1_C, 
                            terms = c("Lat_abs[-1.1:1.7 by=.2]", 
                                      "P_hunt[pursuit]"),
                            ci_level=.85,
                            ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lxh_C2 <-  ggemmeans(TP_mod1_C, 
                            terms = c("Lat_abs[-1.1:2.3 by=.2]", 
                                      "P_hunt[group]"),
                            ci_level=.85,
                            ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lxh_C <- rbind(TP_pred_lxh_C1, TP_pred_lxh_C2)


TPlat_vec_C <- c(max(tempP_modsdf_C2$Lat_abs), mean(tempP_modsdf_C2$Lat_abs),
               min(tempP_modsdf_C2$Lat_abs))

TP_pred_mxl_C <-  ggemmeans(TP_mod2_C, terms = 
                              c("mass_ratio", 
                                "Lat_abs[TPlat_vec_C]"),
                            ci_level=.85,
                            ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(tempP_db$mass_ratio), "scaled:scale"))+ attr(scale(tempP_db$mass_ratio), "scaled:center"),
         lat_real= (as.numeric(as.character(group))* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center")) %>% 
  mutate(across(lat_real, round, 2)) 

Plots

Code

(SP_mass_plot <- ggplot(SP_pred_mass)+
    geom_hline(yintercept = 1, linetype= "dashed", size= 1)+
    geom_ribbon(aes(x= mass_real, y= predicted,
                  ymin= conf.low, ymax= conf.high)
              , alpha=0.7, fill= "#fde725")+
  geom_line(aes(x= mass_real, y=predicted),
            linewidth= 0.8 )+
  labs(x= "ln(Mass ratio)",
        y= "Spatial Overlap",
       tag= "A")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto")))

Code

(SP_phunt_plot <- ggplot(SP_pred_hunt)+
   geom_hline(yintercept = 1, linetype= "dashed", size= 1)+
   geom_errorbar(aes(x= x, y= predicted,
                  ymin= conf.low, ymax= conf.high)
              , size= 1, width= 0.1)+
  geom_point(aes(x= x, y=predicted),
            size= 4)+
  labs(x= "Hunt strategy",
        y= "Spatial Overlap",
       tag= "B")+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"))
  )

Code

# Temporal models

(TP_latxp_plot <- ggplot(TP_pred_lat)+
   geom_ribbon(aes(x= lat_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group,
                  group= group)
              , alpha= 0.7, linewidth= 1)+
    scale_fill_viridis_d()+
  geom_line(aes(x= lat_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "Absolute latitude",
        y= "Temporal Overlap",
       group= "Hunt strategy",
       linetype= "Hunt strategy",
       fill= "Hunt strategy",
       tag= "C")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="vertical",
         legend.position = c(0.7, 0.20),
         legend.background = element_blank()))

Code

(TP_sampxp_plot <- ggplot(TP_pred_sxph)+
   geom_ribbon(aes(x= samp_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group,
                  group= group)
              , alpha= 0.7, linewidth= 1)+
    scale_fill_viridis_d()+
  geom_line(aes(x= samp_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "Survey duration (months)",
        y= "Temporal Overlap",
       group= "Hunt strategy",
       linetype= "Hunt strategy",
       fill= "Hunt strategy",
       tag= "H")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="vertical",
         legend.position = c(0.8, 0.20),
         legend.background = element_blank()))

Code

# Felidae latitude
(TP_lat_plot_F <- ggplot(TP_pred_lat_F)+
   geom_ribbon(aes(x= lat_real, y= predicted,
                  ymin= conf.low, ymax= conf.high),
               alpha=0.7, fill= "#fde725")+
   geom_line(aes(x= lat_real, y=predicted),
            size= 0.8 )+
  labs(x= "Absolute latitude",
        y= "Temporal Overlap",
        tag= "E")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto")))

Code

(TP_mxh_plot_F <- ggplot(TP_pred_mxhunt_F)+
   geom_ribbon(aes(x= mass_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group, 
                  group= group,
                  linetype= group)
              , alpha=0.7)+
    scale_fill_viridis_d()+
  geom_line(aes(x= mass_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "ln(Mass ratio)",
        y= "Temporal Overlap",
       group= "Hunt \nstrategy",
       fill= "Hunt \nstrategy",
       linetype= "Hunt \nstrategy",
       tag= "D")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="horizontal",
         legend.position = c(0.55, 0.10),
         legend.background = element_blank()))

Code

(TP_samp_plot_F <- ggplot(TP_pred_samp_F)+
    geom_ribbon(aes(x= samp_real, y= predicted,
                  ymin= conf.low, ymax= conf.high)
              , alpha=0.7, fill= "#fde725")+
  geom_line(aes(x= samp_real, y=predicted),
            size= 0.8 )+
  labs(x= "Survey duration (months)",
        y= "Temporal Overlap",
       tag= "I")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto")))

Code

(TP_lxh_plot_C <- ggplot(TP_pred_lxh_C)+
   geom_ribbon(aes(x= lat_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group, 
                  group= group,
                  linetype= group)
              , alpha=0.7)+
    scale_fill_viridis_d()+
  geom_line(aes(x= lat_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "Absolute latitude",
        y= "Temporal Overlap",
       group= "Hunt \nstrategy",
       fill= "Hunt \nstrategy",
       linetype= "Hunt \nstrategy",
       tag= "F")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="horizontal",
         legend.position = c(0.55, 0.10),
         legend.background = element_blank()))

Code

(TP_mxl_plot_C <- ggplot(TP_pred_mxl_C)+
   geom_ribbon(aes(x= mass_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= as.factor(lat_real), 
                  group= as.factor(lat_real),
                  linetype= as.factor(lat_real))
              , alpha=0.7)+
    scale_fill_viridis_d()+
  geom_line(aes(x= mass_real, y=predicted, group= as.factor(lat_real),
                linetype= as.factor(lat_real)),
            size= 0.8 )+
  labs(x= "ln(Mass ratio)",
        y= "Temporal Overlap",
       group= "Absolute \nlatitude",
       fill= "Absolute \nlatitude",
       linetype= "Absolute \nlatitude",
       tag= "G")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction= "horizontal",
         legend.position = c(0.55, 0.10),
         legend.background = element_blank()))

Code

(PPrediction_plot <- (SP_mass_plot+ SP_phunt_plot+ TP_latxp_plot+ 
                     TP_mxh_plot_F+ TP_lat_plot_F+ TP_lxh_plot_C+ 
                       TP_mxl_plot_C+ TP_sampxp_plot+ TP_samp_plot_F)+
  plot_layout(ncol = 3))

Code

ggsave(filename = "Figs/Ppreds_plot.png", plot = PPrediction_plot,
        width = 12, height = 10)

References

Arnold, Todd W. 2010. “Uninformative Parameters and Model Selection Using Akaike’s Information Criterion.” The Journal of Wildlife Management 74 (6): 1175–78. https://doi.org/10.2193/2009-367.

Barton, Kamil. 2020. “MuMIn: Multi-Model Inference.” https://CRAN.R-project.org/package=MuMIn.

Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. New York: Springer-Verlag. //www.springer.com/us/book/9780387953649.

Cribari-Neto, Francisco, and Achim Zeileis. 2010. “Beta Regression in r” 34. https://doi.org/10.18637/jss.v034.i02.

Faurby, Søren, Rasmus Ø Pedersen, Matt Davis, Simon D. Schowanek, Scott Jarvie, Alexandre Antonelli, and Jens-Christian Svenning. 2020. MegaPast2Future/PHYLACINE_1.2: PHYLACINE Version 1.2.1. Zenodo. https://doi.org/10.5281/zenodo.3690867.

Hartig, Florian. 2021. “DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models.” https://CRAN.R-project.org/package=DHARMa.

Hassanin, Alexandre, Géraldine Veron, Anne Ropiquet, Bettine Jansen van Vuuren, Alexis Lécu, Steven M. Goodman, Jibran Haider, and Trung Thanh Nguyen. 2021. “Evolutionary History of Carnivora (Mammalia, Laurasiatheria) Inferred from Mitochondrial Genomes.” PLOS ONE 16 (2): e0240770. https://doi.org/10.1371/journal.pone.0240770.

Lüdecke, Daniel, Mattan S. Ben-Shachar, Indrajeet Patil, Philip Waggoner, and Dominique Makowski. 2021. “Performance: An r Package for Assessment, Comparison and Testing of Statistical Models” 6: 3139. https://doi.org/10.21105/joss.03139.

Makowski, Dominique, Mattan S. Ben-Shachar, Indrajeet Patil, and Daniel Lüdecke. 2020. “Estimation of Model-Based Predictions, Contrasts and Means.” https://github.com/easystats/modelbased.

Mazerolle, Marc J. 2023. “AICcmodavg: Model Selection and Multimodel Inference Based on (q)AIC(c).” https://cran.r-project.org/package=AICcmodavg.

Revelle, William. 2022. “Psych: Procedures for Psychological, Psychometric, and Personality Research.” https://CRAN.R-project.org/package=psych.

Richmond, Orien M. W., James E. Hines, and Steven R. Beissinger. 2010. “Two-Species Occupancy Models: A New Parameterization Applied to Co-Occurrence of Secretive Rails.” Ecological Applications 20 (7): 2036–46. https://doi.org/10.1890/09-0470.1.

Ridout, M. S., and M. Linkie. 2009. “Estimating Overlap of Daily Activity Patterns from Camera Trap Data.” Journal of Agricultural, Biological, and Environmental Statistics 14 (3): 322–37. https://doi.org/10.1198/jabes.2009.08038.

Schloerke, Barret, Di Cook, Joseph Larmarange, Francois Briatte, Moritz Marbach, Edwin Thoen, Amos Elberg, and Jason Crowley. 2021. “GGally: Extension to ’Ggplot2’.” https://CRAN.R-project.org/package=GGally.

Soria, Carmen D., Michela Pacifici, Moreno Di Marco, Sarah M. Stephen, and Carlo Rondinini. 2021. “COMBINE: A Coalesced Mammal Database of Intrinsic and Extrinsic Traits.” Ecology 102 (6): e03344. https://doi.org/10.1002/ecy.3344.

Sutherland, C., D. Hare, P. J. Johnson, D. W. Linden, R. A. Montgomery, and E. Droge. 2023. “Practical Advice on Variable Selection and Reporting Using AIC.” Proceedings of the Royal Society B: Biological Sciences. https://ora.ox.ac.uk/objects/uuid:38650aad-20a6-4980-af4f-7ded0968e908?fbclid=IwAR29QmTXxKFknl9BI_SKw-6I5gL-9NmdaytggzFomgMOPMY9V50q9ohJQ6k.

Waddle, J.Hardin, Robert M. Dorazio, Susan C. Walls, Kenneth G. Rice, Jeff Beauchamp, Melinda J. Shuman, and Frank Mazzotti. 2010. “A New Parameterization for Estimating Co-Occurrence of Interacting Species.” Ecological Applications 20 (5): 302315. https://doi.org/https://doi.org/10.1890/09-0850.1.

Zuur, Alain F., Elena N. Ieno, and Chris S. Elphick. 2010. “A Protocol for Data Exploration to Avoid Common Statistical Problems.” Methods in Ecology and Evolution 1 (1): 3–14. https://doi.org/https://doi.org/10.1111/j.2041-210X.2009.00001.x.

Zuur, Alain F, Elena N Ieno, Neil J Walker, Anatoly A Saveliev, and Graham M Smith. 2009a. Mixed Effects Models and Extensions in Ecology with r. Vol. 574. Springer.

———. 2009b. Mixed Effects Models and Extensions in Ecology with r. 1st ed. Statistics for Biology and Health. Springer.